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Abstract. We present a detailed analysis of the modulational instability of the 
zone-boundary mode for one and higher-dimensional Fermi-Pasta-Ulam (FPU) 
lattices. The growth of the instability is followed by a process of relaxation to 
equipartition, which we have called the Anti-FPU problem because the energy is 
initially fed into the highest frequency part of the spectrum, while in the origi¬ 
nal FPU problem low frequency excitations of the lattice were considered. This 
relaxation process leads to the formation of chaotic breathers in both one and 
two space dimensions. The system then relaxes to energy equipartition, on time 
scales that increase as the energy density is decreased. We supplement this study 
by considering the nonconservative case, where the FPU lattice is homogeneously 
driven at high frequencies. Standing and travelling nonlinear waves and solitonic 
patterns are detected in this case. Finally we investigate the dynamics of the FPU 
chain when one end is driven at a frequency located above the zone boundary. We 
show that this excitation stimulates nonlinear bandgap transmission effects. 


1 Introduction 


In 1955, reporting about one of the first numerical simulations, Fermi, Pasta and Ulam (FPU) [T] 
remarked that it was ... very hard to observe the rate of “thermalization” or mixing ... in a 
nonlinear one-dimensional lattice in which the energy was initially fed into the lowest frequency 
mode. Even if the understanding of this problem has advanced significantly afterwards na, 
several issues are still far from being clarified. In most cases, the evolution towards energy 
equipartition among linear modes has been checked considering an initial condition where all 
the energy of the system was concentrated in a small packet of modes centered around some 
low frequency. 

Beginning with the pioneering paper of Zabusky and Deem [1], the opposite case in which 
the energy is put into a high frequency mode has been also analyzed. In this early paper, the 
zone-boundary mode was excited with an added spatial modulation for the one-dimensional 
a-FPU model (quadratic nonlinearity in the equations of motion). Here, we will study the time- 
evolution of this mode without any spatial modulation for the /5-FPU model (cubic nonlinearity 
in the equations of motion) and some higher-order nonlinearities. Moreover, we will extend the 
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study to higher dimensional lattices. Since the energy is fed into the opposite side of the linear 
spectrum, we call this problem the Anti-FPU problem. 

In a paper by Bundinsky and Bountis [5], the zone-boundary mode solution of the one¬ 
dimensional FPU lattice was found to be unstable above an energy threshold which scales like 
1/A^, where N is the number of oscillators. This result was later and independently confirmed 
by Flach [6] and Poggi et al [7], who also obtained the correct factor in the large iV-limit. These 
results were obtained by a direct linear stability analysis around the periodic orbit corresponding 
to the zone-boundary mode. Similar methods have been recently applied to other modes and 
other FPU-like potentials by Chechin et al [819] and Rink [10]. 

A formula for Ec, valid for all N, has been obtained in Refs. |1:II12I11I14| in the rotating 
wave approximation, and will be also discussed in this paper. Associated with this instability 
is the calculation of the growth rates of mode amplitudes. The appropriate approach for Klein- 
Gordon lattices was first introduced by Kivshar and Peyrard |15j . following an analogy with 
the Benjamin-Feir instability in fluid mechanics m- 

Previously, a completly different approach to describe this instability was introduced by 
Zakharov and Shabat [HI, studying the associated Nonlinear Schrodinger equation in the con¬ 
tinuum limit. A value for the energy threshold was obtained in Ref. [18] in the continuum limit. 
The full derivation starting from the FPU equation of motions was then independently obtained 
by Berman and Kolovskii m in the so-called “narrow-packet” approximation. 

Only very recently the study of what happens after the modulational instability develops has 
been performed for Klein-Gordon [20] and FPU-lattices |21ll3] . From these analyses it turned 
out that these high-frequency initial conditions lead to a completely new dynamical behavior in 
the transient time preceeding the final energy equipartition. In particular, the main discovery 
has been the presence on the lattice of sharp localized modes [HEo]. These latter papers were 
the first to make the connection between energy relaxation and intrinsic localized modes [22] . 
or breathers m- Later on, a careful numerical and theoretical study of the dynamics of a 
/3-FPU model was performed [H] . It has been shown that moving breathers play a relevant role 
in the transient dynamics and that, contrary to exact breathers, which are periodic solutions, 
these have a chaotic evolution. This is why they have been called chaotic breathers. Following 
these studies, Lepri and Kosevich [^ and Lichtenberg and coworkers have further 

characterized the scaling laws of relaxation times using continuum limit equations. 

Zabusky et al. have recently simulated numerically the behavior of the one-dimensional, 
periodic a-FPU model with optical and acoustic initial excitations of small-but finite and large 
amplitudes. Using beautiful color representations [28] of the numerical results, they find nearly 
recurrent solutions, where the optical result is due to the appearance of localized breather-like 
packets. For large amplitudes, they obtained also chaotic behaviors for the alpha lattice. 

Using a theory (originally developed [29180] for the discrete nonlinear Schrodinger equation) 
where standard Gibbsian equilibrium statistical mechanics was considered to predict macro¬ 
scopic average values for a thermalized state in the thermodynamic limit, Johansson has re¬ 
cently analyzed m certain aspects of a mixed Klein-Gordon/FPU chain. In particular, he 
shows that the available phase space is divided into two separated parts with qualitatively dif¬ 
ferent properties in thermal equilibrium; one part corresponding to a normal thermalized state 
with exponentially small probabilities for large-amplitude excitations, and another part typ¬ 
ically associated with the formation of high-amplitude localized excitations, interacting with 
an infinite-temperature phonon bath. Observing the /3-FPU chains in the thermalized state, 
Gershgorin et al showed [32] via numerical simulation that discrete breathers actually persist 
and have a turbulent-like behavior. They describe the dynamical scenario as spatially highly 
localized discrete breathers riding chaotically on spatially extended, renormalized waves. 

Recently Flach et al [53] have focused on the main FPU observation that the initially excited 
normal mode shares its energy for long times only with a few other modes from a frequency 
neighbourhood in modal space. They have identified this long lasting regime as a dynamical 
localization effect and applied the methods developed for discrete breathers in FPU chains 
to the dynamics of normal modes. The result is that time-periodic and modal-space-localized 
orbits, (they call q-breathers) persist in the FPU model. The dynamics generated by one initially 
excited mode evolves close to the related q-breathers for very long times. Thus many features 
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of the short- and medium-time evolution of natural packets are encoded in the prohle of these 
objects. 

Let us briefly comment about the chaoticity of these spontaneously created breathers. 
Tailleur and Kurchan have recently implemented m an efficient method that allows one to se¬ 
lect trajectories with unusual chaoticity, with Lyapunov weighted dynamics (LWD) (a method 
originally proposed in the context of chemical reactions |35jl. As an example of application, 
they study the Fermi-Pasta-Ulam nonlinear chain starting from a microcanonical equilibrium 
configuration. They show that the algorithm rapidly singles out the chaotic-breathers when 
searching for trajectories with high level of chaoticity (typically they study cases where the 
Lyapunov is three times the one of a typical equilibrium trajectory), thus confirming that the 
large Lyapunov conhgurations are dominated by chaotic breathers. 

Most of the previous studies are for one-dimensional lattices. We have recently derived mod- 
ulational instability thresholds also for higher dimensional lattices |36j and we have presented 
a study of chaotic breathers formation in two-dimensional FPU lattices. However, the study is 
extremely preliminary and further analyses are needed. In particular, the full process of relax¬ 
ation to energy equipartion and the associated time scales have not been carefully studied in 
two-dimensional FPU lattices. Pioneering results on the relaxation process in a two-dimensional 
triangular FPU lattice from low frequency initial states seem to indicate a faster evolution to 
equipartition m- Benettin has discovered that for large values of the energy per site the time 
scale for equipartition can be quite short, even in the thermodynamic limit of the lattice size, 
and that this time scale increases as only one-over the energy per site. If one lowers the energy 
per site below a critical threshold, however, the time scale for equipartition on finite lattices 
grows much more rapidly. But the critical threshold value of energy per site appears to vanish as 
the lattice size goes to infinity. A similar analysis for high frequencies remains to be performed. 

Further, having already studied the process of formation of stable localized structures arising 
from modulational instability in the conservative case EH, we are strongly motivated to see how 
the presence of forcing and damping affects this process. To remain close to the Hamiltonian 
case, we restrict ourselves to the case of small damping. Various types of forcing are in principle 
possible, depending on the physical situation under study. However, a general requirement for 
localization is to excite band-edge modes. For Klein-Gordon lattices this is naturally realized 
using a spatially uniform driving field, which has been shown to induce interesting pattern for¬ 
mation phenomena m- On the other hand, this forcing would not be effective for FPU lattices, 
because, due to the symmetry of the Hamiltonian, the zero mode is decoupled. Alternatively, 
since spatial localization appears from the instability of band-edge modes, we choose in the 
context of Anti FPU scenario to drive the system near the zone boundary wavelength. As it 
will be shown below stationary localized patterns (either moving or static) appear under such 
a homogeneous driving and damping. 

Finally, we consider driving the system by one end, again with frequencies above the zone 
boundary. By this we stimulate the appearance of a “supratransmission” scenario [38] in FPU, 
i.e. the chain becomes conductive only for driving amplitudes above a threshold [35]. This 
phenomenon is explained in terms of nonlinear response manifolds in Ref. |40j . 

We have organized the paper in the following way. In Sectional the modulational instability 
of zone-boundary modes on the lattice is discussed, beginning with the one-dimensional case, 
followed by the two-dimensional and higher dimensional cases and finishing with the continuum 
nonlinear Schrbdinger approach. In this Section, we also decribe the mechanism of creation of 
chaotic breathers in one and two dimensions. Section|3|deals with the driven-damped Anti-FPU 
scenario (homogeneous and point-like driving). Some final remarks and conclusions are reported 
in Section m 


2 Modulational Instability 

2.1 The one-dimensional case 

We will discuss in this section modulational instability for the one-dimensional FPU lattice, 
where the linear coupling is corrected by a (2i + l)th order nonlinearity, with I a positive 
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integer. Denoting by the relative displacement of the n-th particle from its equilibrium 

position, the equations of motion are 

tin — t/n+1 “t” tin —r “b (ttn+t tin) (ttn tin —r) • (1) 

We adopt a lattice of N particles and we choose periodic boundary conditioirs. For the sake 
of simplicity, we first report on the analysis for £ = 1 (i.e. for the /3-FPU model) and then we 
generalize the results to any £-value. 

Due to periodic boundary conditions, the normal modes associated to the linear part of 
Eq. ([I} are plaire waves of the form 

= (2) 

where On{t) = qn — Lot and q = 2TYk/N {k = —N/2,...,N/2). The dispersion relation of 
nonlinear phonons in the rotating wave approximation m is io'^iq) = 4(1 + a) sin^(g/2), where 
a = 3a^sin^(q/2) takes into account the nonlinearity. Modulational instability of such a plane 
wave is investigated by studying the linearized equation associated with the envelope of the 
carrier wave ©. Therefore, one introduces infinitesimal perturbations in the amplitude and 
phase and looks for solutions of the form 

Unit) = ^[1 + b^it)] + ^[1 + 6„(t)] 

= a[l + 6„(f)] cos[qn - Lot + (3) 

where and il^n are reals and assumed to be small in comparison with the parameters of 
the carrier wave. Substituting Eq. m into the equations of motion and keeping the second 
derivative, we obtain for the real and imaginary part of the secular term the following 

equations 

- Lo'^bn + 2i0ipn + = (1 + 2a) [cos q (6„+i + 6„_i) - 26„] 

- a (6„+i + bn-i - 2bnCOSq) - (1 + 2a) sinq (ipn+i - i^n-i) (4) 
-Lo’^i^n - 2Lobn + '0„ = (1 + 2a) [cos q (V'n +1 + t/'n-l) “ 2^„] 

+ (1 + 2a) sing (6„+i - 6„_i) + a (t/’n+i + ipn-i - 2ipn cosq). (5) 

Further assuming = bo c.c. and ipn = "00 _|_ q q obtain the 

following equations for the secular term 


+co^ + 2(1 + 2a) (cos g cos Q — 1) — 2a(cosQ — cosg) 


— 2iil)o [ujfl + (1 + 2a) sing sin Q] = 0 


l/’o 


+Lo^ + 2(1 + 2a)(cosgcos(5 — 1) + 2a(cosQ — cosg) 


+ 2ibo [iofi + (1 + 2a) sin g sin Q] = 0. 


( 6 ) 

(7) 


In the case of Klein-Gordon type equations |15I20) . one neglects the second order derivatives in 
Eqs. dH)-®. This can be justified by the existence of a gap in the dispersion relation for g = 0, 
which allows to neglect 17^ with respect to In the FPU case, this approximation is worse, 
especially for long wavelengths, because there is no gap. 

Non trivial solutions for Eqs. can be found only if the Cramer’s determinant vanishes, 

i.e. if the following equation is fulfilled: 


(17 + w)^ — 4(1 + 2a) sin^ 


g + <5 


(17 — w)^ — 4(1 + 2a) sin^ 


= 4a^ (cos Q — cos g)^ 


q-Q 


(8) 
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This equation admits four different solutions when the wavevectors q of the unperturbed 
wave and Q of the perturbation are fixed. If one of the solutions is complex, an instability of one 
of the modes {q ± Q) is present, with a growth rate equal to the imaginary part of the solution. 
Using this method, one can derive the instability threshold amplitude for any wavenumber. A 
trivial example is the case of g = 0, for which we obtain fl = ± sin ((5/2), which proves that 
the zero mode solution is stable. This mode is present due to the invariance of the equations 
of motion (HD with respect to the translation + const and, as expected, is completely 

decoupled from the others. 

A first interesting case is q = tt. One can easily see that Eq. (|SD admits two real and two 
complex conjugate imaginary solutions if and only if 


cos^ 



1 + a 
1 + 3a 


(9) 


This formula was first obtained by Sandusky and Page (Eq. (22) in Ref. [T3]) using the rotating 
wave approximation. The first mode to become unstable when increasing the amplitude a 
corresponds to the wavenumber Q = 2tt/N. Therefore, the critical amplitude Oc above which 
the q = TT-mode looses stability is 


/ sin^ \ 

\3 [3 cos^ (tt/IV) — 1] y 

This formula is valid for all even values of N and its large A^-limit is 


( 10 ) 


In Fig. [U we show its extremely good agreement with the critical amplitude determined from 
numerical simulations. It is interesting to emphasize that the analytical formula (HHD diverges 
for N = 2, predicting that the 7r-mode is stable for all amplitudes in this smallest lattice. This 
is in agreement with the Mathieu equation analysis (see Ref. [7] p. 265). 

It is also interesting to express this result in terms of the total energy to compare with what 
has been obtained using other methods |17| 5 |19| 6 |7j. Since for the 7r-mode the energy is given 
by E = N{2a? + da"'), we obtain the critical energy 


2N . 2 ( '’^ \ 7 cos^ (^r/fV) — 1 
9 \N) [3 cos^ (tt/IV) — 1]^ 


( 12 ) 


For large N, we get 

^ ^ (t^) ■ 

This asymptotic behavior is the same as the one obtained using the narrow packet approxima¬ 
tion in the context of the nonlinear Schrddinger equation by Berman and Kolovskii (Eq. (4.1) 
in Ref. m)- The correct scaling behavior with N of the critical energy has been also obtained 
by Bundinsky and Bountis (Eq. (2.22) in Ref. [5]) by a direct linear stability analysis of the tt- 
mode. The correct formula, using this latter method, has been independently obtained by Flach 
(Eq. (3.20)) in Ref. [5]) and Poggi and Ruffo (p. 267 of Ref. [7]). Recently, the iV“'-scaling of 
formula (USD has been confirmed using a different numerical method and, interestingly, it holds 
also for the 2tt/3 and 7r/2 modes |41) . 

This critical energy is also very close to the Chirikov “stochasticity threshold” energy ob¬ 
tained by the resonance overlap criterion for the zone boundary mode [42]. The stochasticity 
threshold phenomenon has been thoroughly studied for long wavelength initial conditions, and 
it has been clarified that it corresponds to a change in the scaling law of the largest Lyapunov 
exponent [43] . We will show in Section [3] that above the modulational instability critical energy 
for the TT-mode one reaches asymptotically a chaotic state with a positive Lyapunov exponent, 
consistently with Chirikov’s result. 
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The above results can be generalized to nonlinearities of + 1 order in the equations of 
motion O- We limit the analysis to the rr-mode, for which the instability condition Q takes 
the form 


2 Q ^ 1 + a 

2 l + (2£+l)a’ 


(14) 


where 


a = 


( 2 ^ + 1 )! 2e 
£!(£ + !)! 


(15) 


Hence the critical amplitude above which the rr-mode is unstable is 


dc 


£!(£+!)! sin^ {n/N) 
(2£+l)![(2£+l)cos2 (tt/N) 



(16) 


leading to the large N scaling 

(17) 

(18) 


a, - 
Ee ~ 


This scaling also corresponds to the one found in Ref. [44] when discussing tangent bifurcations 
of band edge plane waves in relation with energy thresholds for discrete breathers. Their “de¬ 
tuning exponent” z has a direct connection with the nonlinearity exponent £ = zj^. We will 
see in Section that this analogy extends also to higher dimensions. 

For fixed N, Qc is an increasing function of the power of the coupling potential with the 
asymptotic limit lim^_^oo Oc = 0.5. Therefore, in the hard potential limit the critical energy 
for the TT-mode increases proportionally to N. The fact that we find a higher energy region 
where the system is chaotic is not in contradiction with the integrability of the one-dimensional 
system of hard rods [4^, because in the present case we have also a harmonic contribution at 
small distances. 

For the FPU-a model (quadratic nonlinearity in the equations of motion), the 7r-mode is 
also an exact solution which becomes unstable at some critical amplitude which, contrary to 
the case of the FPU-/3 model, is Windependent [T3I8) : which means that the critical energy is 
proportional to N and then that 7r-mode can be stable in some low energy density limit also in 
the thermodynamic limit. 

It has also been realized |7l47l46l8l48lin] that group of modes form sets which are invariant 
under the dynamics. The stability analysis of pair of modes has shown a complex depen¬ 

dence on their relative amplitudes. The existence of such invariant manifolds has also allowed 
to construct Birkhoff-Gustavson normal forms for the FPU model, paving the way to KAM 
theory [49]. 


2.2 Higher dimensions 

In this Section, we will first discuss modulational instability of the two-dimensional FPU model. 
The method presented in Section 12.11 can be easily extended and the global physical scenario 
is preserved. However, the scaling with N of the critical amplitude changes in such a way to 
make critical energy constant, in agreement with the analysis of Ref. [44] . 

The masses lie on a two-dimensional square lattice with unitary spacing in the (x, y) plane. 
We consider a small relative displacement Un,m {n,m € [IjA/^]) in the vertical direction z. 
Already with an harmonic potential, if the spring length at equilibrium is not unitary, the 
series expansion in Un,m of the potential contains all even powers. We retain only the first 
two terms of this series expansion. After an appropriate rescaling of time and displacements 
to eliminate mass and spring constant values, one gets the following adimensional equations of 
motions 

fln,m — 4“ ^n —l,m 4“ 4- Un,m — 1 

4- (rin+l.m 4- (r^n—l,m '^n,m) 4- {Un^m+1 4- (rin,m—1 '^n^m )(19) 
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Fig. 1. Left panel: Modulational instability threshold amplitude for the rr-mode versus the number of 
particles in the one-dimensional FPU lattice. The solid line corresponds to the analytical formula m, 
the dashed line to its large Westimate (El) and the diamonds are obtained from numerical simula¬ 
tions. Right panel: Modulation instability threshold for the (tt, tt) mode versus number of oscillators in 
two dimensional array. Solid line is derived from the exact analytical consideration ([5U1). dashed line 
describes the estimate from NLS equation in large N limit (O and diamonds are result of numerical 
simulations. 


Considering periodic boundary conditions, plane waves solutions have the form 


Un,m = a cos {qxn + qym — wt). 

In the rotating wave approximation[T3], one immediately obtains the dispersion relation 


2 A ■ 2 , A ■ 2 % 

o;^ = 4 sm-h 4 sm — 

2 2 


12a^ 


sm 


4 Qx 


sm 


4 ^ 

2 J 


( 20 ) 

( 21 ) 


which becomes exact for the zone-boundary mode (q^^qy) = (tt,tt), 

<, = 8(l + 3a2). (22) 

To study the stability of the zone-boundary mode, we adopt a slightly different approach. 
Namely, we consider the perturbed relative displacement field of the form 

Un,m = (^ + + c.C., (23) 

where is complex. This approach turns out to be equivalent to the one of Section 12.11 in 
the linear limit. 

Substituting this perturbed displacement field in Eqs. m, we obtain 


[1 + 2a] [bn+i ,m ^n,m — 1 

[bn+l^m + ^n-l.m + ^n,m-|-l + „j] = —bn,m + ‘2‘iuJ-x^-nbn^m + W-n-,'7r^n,rr(?4) 

where a = 3a^. Looking for solutions of the form 

^ _ j^^i{Qxn+Qym-Qt) _|_ j^^-i(Qxn+Qym-f2t) ^"2^^ 

we arrive at the following set of linear algebraic equations for the complex constants A and B 

[(C -h a;^,^)2 - 8(1 -h 2a)A] A + 8aAB = 0 (26) 

8aAA + [(17 - - 8(1 -h 2a)Z\] B = 0, (27) 
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where 2A = cos^{Qx/‘2) + cos^{Qy/2). As for the one-dimensional case, we require that the 
determinant of this linear system in A and B vanishes, which leads to the following condition 

[{n + - 8Z\(1 -h 2a)] [(f? - - 82\(1 -h 2a)] = bda^Zl^. (28) 

This equation admits two real and two complex conjugated imaginary solutions in 12 if 

1 -I- a 


Z\ > 


l-b3a 


(29) 


which is the analogous of condition (jH]) for two dimensions. One can achieve the minimal nonzero 
value of the r.h.s. of the above expression choosing Qx = 0, Qy = 2tt/N, which leads to the 
following result for the critical amplitude 


/ sin^(7r/A^) \ 

\3[3cos^(7r/iV) -1-1]/ 


(30) 


Its large N limit is 


(31) 


This prediction is compared with numerical data in Fig. [T] The agreement is good for all values 
of N. 

Since the relation between energy and amplitude is now E = 2iV^(2a^ -|-4a^), we obtain the 
critical energy in the large A^-limit as 



(32) 


This shows that the critical energy is now constant in the thermodynamic limit, which agrees 
with the remark of Ref. [H] about the existence of a minimal energy for breathers formation . 

The results of this Section can be easily extended to any dimension d. Repeating the same 
argument, we arrive at the following estimates for the critical amplitude and energy in the large 
N limit 

Ea = + O . (34) 

O 

This means that the critical energy density £c = E^/N for destabilizing the zone boundary 
mode vanishes as 1/A^^, independently of dimension. 


2.3 Large N limit using the Nonlinear Schrodinger equation 


The large N limit expressions (13311 and (I34|) can be derived also by continuum limit considera¬ 
tions. We will derive the general expression for any dimension d. The displacement field can be 
factorized into a complex envelope part multiplied by the zone boundary mode pattern in d 
dimensions. 






. 


(35) 


where 

w^,...,7r = Vd[4(l -I- 3|^/p)]. (36) 

Substituting Eq. dsni) into the FPU lattice equations in d dimensions, a standard proce¬ 
dure |51I50) leads to the following d dimensional Nonlinear Schrodinger (NTS) equation: 



+ ^ - QV’1'01^ = 0, 


(37) 
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where Ad is the d dimensional Laplacian. The parameters P and Q are derived from the 
nonlinear dispersion relation 


i=l 


[4sin2| + 12|V^|2sin4:^ 


2 J 


(38) 


as 


P = (gi = TT, . . . , gd = TT, \tlj\ = 0) = 


d'^uj 

duj 


2y/d 


Q = (gi = TT,..., gd = TT, \ip\ = 0) = -3^. 


(39) 

(40) 


Assuming that, at the first stage, modulation instability develops along a single direction 
X and that the field remains constant along all other directions, one gets the one-dimensional 
NTS equation 


.d'ip ^ P 


Qm^ = o. 


(41) 


Following the results of the inverse scattering approach m, any initial distribution of amplitude 
lipl and length A along x, and constant along all other directions, produces a final localized 
distribution if m 




p 

Q 


(42) 


This means that if the initial state is taken with constant amplitude |'0| = a on the d-dimensional 
lattice with N'^ oscillators, the modulation instability threshold is 


{acNf 


6d 


(43) 


which coincides with the leading order inEq. dMl). 


3 Emergence of Localizations in Anti-FPU 

3.1 Conservative Case: Chaotic Breathers 

In this Section, we will discuss what happens when the modulational energy threshold is over¬ 
come. The first thorough study of this problem can be found in Ref. many years after the 
early pioneering work of Zabusky and Deem [4] . Already in Ref. [21] , it has been remarked that 
an energy localization process takes place, which leads to the formation of breathers [23]. This 
process has been further characterized in terms of time-scales to reach energy equipartition and 
quantitative localization properties in Ref. |24| . The localized structure which emerges after 
modulational instability has been here called “chaotic breather” (CB). The connection between 
CB formation and continuum equations has been discussed in Refs. ESllZ], while the relation 
with the process of relaxation to energy equipartition has been further studied in Ref. |26| . We 
will briefly recall some features of the localization process in one dimension and present new 
results for two dimensions. 

For long time simulations, we use appropriate symplectic integration schemes in order to 
preserve as far as possible the Hamiltonian structure. For the one dimensional FPU, we adopt 
a 6th-order Yoshida’s algorithm m with a time step dt = 0.01; this choice allows us to 
obtain an energy conservation with a relative accuracy AE/E ranging from 10“^° to 10“^^. 
For two dimensions, we use instead the 5-th order symplectic Runge-Kutta-Nystrom algorithm 
of Ref. |53] , which gives a similar quality of energy conservation. 
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t3=10® (d) 



t2=3.10‘‘ (c) 



n 


tj = 500 (t>) 



Fig. 2. Time evolution of the local energy (gH). In panel (a), the horizontal axis indicates lattice sites 
and the vertical axis is time. The grey scale goes from En = 0 (white) to the maximum i?„-value 
(black). The lower rectangle corresponds to 0 < t < 3000 and the upper one to 5.994 10® < t < 6.10®. 
Figs, (b), (c) and (d) show the instantaneous local energy En along the N = 128 chain at three different 
times. Remark the difference in vertical amplitude in panel (c), when the CB is present. The initial 
TT-mode amplitude is a = 0.126 > — 0.010. 


We report in Fig. UKa) a generic evolution of the one dimensional 7r-mode above the mod¬ 
ulation instability critical amplitude (a > Oc). The grey scale refers to the energy residing on 
site n, 

'^n) T 2^^^^ tXyj_i), (44) 

where the FPU-potential is V{x) = -I- \x'^. Figs. [2](b), [2Kc) and [2](d) are three successive 

snapshots of the local energy En along the chain. At short time, a slight modulation of the 
energy in the system appears (see Fig. [2Kb)) and the 7r-mode is destabilized [13]. Later on, 
as shown in Fig. UKa), only a few localized energy packets emerge: they are breathers [23] . 
Inelastic collisions of breathers have a systematic tendency to favour the growth of the big 
breathers at the expense of small ones [54155] . Hence, in the course of time, the breather 
number decreases and only one, of very large amplitude, survives (see Fig. [2Kc)): this is the 
localized excitation we have called chaotic breather (CB). The CB moves along the lattice with 
an almost ballistic motion: sometimes it stops or reflects. During its motion the CB collects 
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Fig. 3. Panel (a) presents the evolution of Co(t) of formula (1451) for the one-dimensional FPU lattice 
with N = 128 oscillators, initialized on the 7r-mode with an amplitude a — 0.126 > — 0.010. The 

dashed line indicates the equilibrium value Co — 1.795. Panel (b) presents the corresponding finite time 
largest Lyapunov exponent. Panel (c) shows Co{t) for the two-dimensional FPU lattice with 20 20 

oscillators, initialized on the (7r,7r)-mode with an amplitude a = 0.425 > — 0.045. Panel (d) presents 

the finite time largest Lyapunov exponent for two dimensions. 


energy and its amplitude increases. It is important to note that the CB is never at rest and that 
it propagates with a given subsonic speed m- Finally, the CB decays and the system reaches 
energy equipartition, as illustrated in Fig. [2jd). 

In order to obtain a quantitative characterization of energy localization, we introduce the 
“participation ratio” 

N 

Co{t)=N (45) 

(S'l 

which is of order one if Ei = E/N at each site of the chain and of order N if the energy 
is localized on only one site. In Fig. EKa), Cq is reported as a function of time. Initially, Cq 
grows, indicating that the energy, evenly distributed on the lattice at t = 0, localizes over a 
few sites. This localized state survives for some time. At later times, Cq starts to decrease and 
finally reaches an asymptotic value Cq which is associated with the disappearance of the CB 
(an estimate of Cq has been derived in Ref. |24j taking into account energy fluctuations and is 
reported with a dashed line in Fig.[2Ka)). At this stage, the energy distribution in Fourier space 
is flat, i.e. a state of energy equipartition is reached. 

In Fig. ISKb), we show the finite time largest Lyapunov exponent Xi{t) for the same orbit 
as in Fig. EKa). We observe a growth of Xi{t) when the CB emerges on the lattice and a 
decrease when it begins to dissolve. The peak in Ai(t) perfectly coincides with the one in Cq- 
Due to this increase of chaos associated with localization, we have called the breather chaotic 
(although chaos increase could be the result of more complicated processes of interaction with 
the background). 
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In Ref. [21] , the time-scale for the relaxation to equipartition has been found to increase as 
{E/N) in the small energy limit. This has been confirmed by the followers of this study [25126127] . 

Such power law scalings are found also for the FPU relaxation starting from long wave¬ 
lengths m-- the so-called FPU problem. We have termed the relaxation process which starts 
from short wavelengths the Anti-FPU problem., just because of the similarities in the scaling 
laws. The main feature of the latter problem is that relaxation to equipartition goes through 
a complex process of localized structures formation well described by breathers or, in the low- 
amplitude limit, by solitons of the NonLinear Schrodinger equation. On the contrary, for the 
original FPU problem, an initial long wavelength excitation breaks up into a train of mKdV- 
solitons. The final relaxation to equipartition is however due to an energy diffusion process 
which has similar features for both the FPU and the anti-FPU problem [2S]. 

A similar evolution of the local energy 


E — 

^n,m — 


T^U('Un+l,m ’^n,rri} T ^ l^('an,m-t-l ^n,m) 

— ‘^n,m) T ^f^('^7i,77i—1 ‘^n,rn) 


(46) 


is observed for the two-dimensional case (see Fig.|4|). In this figure, we just show the initial evo¬ 
lution which leads to the breathers formation. As for the one-dimensional case, bigger breathers 
eat up smaller ones, and finally only two breathers survive. We don’t observe coalescence to 
a single breather because collisions are more rare in two dimensions. After the formation of 
a few localized structures, one also observes the final relaxation to equipartition which is not 
shown in Fig.|4| This latter is instead evident from the time evolution of C'o(t), the localization 
parameter, shown in Fig. me) : its behavior is very similar to the one-dimensional case. Indeed, 
also the largest finite time Lyapunov exponent behaves similarly (see Fig. EKd)). 


t = 28.5 


t = 533.7 



Fig. 4. Local energy (M surface plots for the two-dimensional FPU lattice with 20 * 20 oscillators, 
initialized on the (tt, 7r)-mode with an amplitude a = 0.225496 > Uc — 0.0453450. Snapshots at four 
different times t are shown. Breathers form after a coalescence process similarly to the one-dimensional 
case. The mobility of the breathers is evident and one also observes in the last panel the final decrease. 
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lattice site 


Fig. 5. Left panel: Static multibreather pattern generated after the modulational instability for p — tt 
{lu = 2.4, / = 0.2250). This pattern stabilizes at t ~ 510^. Right panel: Travelling multibreather 
pattern of the p ^ vr-mode for a lattice of = 512 sites. Here, lo = 2.05, / = 0.075 and p = 2.4542. 


3.2 Nonconservative Case: Homogeneously Driven-Damped Anti-FPU 


The equations of motion of the “externally driven” and damped FPU chain read as follows: 

Un = Un+1 + Un-1 - 2Un + (Wn+l - Unf d" (Wn-l “ + / COs(a;t -f pu), (47) 


where the forcing and damping strengths are gauged by the parameters / and 7 , respectively; 
u} and p are the driving frequency and wavenumber. Considering the Anti-FPU situation, we re¬ 
strict ourselves to the case 7r/2 < |p| < tt. Moreover, here we present only the results concerning 
the range of driving frequencies |a;| > ujp = \/2{l — cosp) for which stationary multibreather 
states develop (they are static if p = tt and move for other cases). For other driving frequencies, 
one deals with the stationary periodic patterns described in Refs. |64I65] . 

Examples of appearance of either static (p = tt) or traveling (p ^ tt) multibreather states 
are shown in Fig.[5j where we plot the local energy vs. the lattice position sampled at the period 
of the forcing. The corresponding spatial Fourier spectrum is shown in Fig. [51 The broad band 
structure of the spectrum reflects the non perfect periodic arrangement of the localized peaks 
in Fig. [3 

Such states can be described in terms of soliton solutions of an associated suitable driven- 
damped nonlinear Schrodinger (NTS) equation. Let us first make the following definition; 


ap{n, -b a+(n, 


(48) 


where Op and its conjugate are smooth functions of n. Such an assumption is possible if the 
wavepacket is concentrated around the driving mode Ak <g; tt. Substituting (l4^ into the 
equations of motion (H71) . one gets: 


da^ 


da 


(w^ - u;'^)ap +2iuj ( ^ 


aip d^ap 3 4 I |2 . j. 


(49) 


t' = 


2 2 


2u} 




uj'^ — 


LU-n 


[n — vt), 


After performing the following re-scalings 
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4.0 

Ek 


2.0 


0.0 



Fig. 6. Left panel: Spatial spectrum of the travelling multibreather, ek = \Uk\ + uJ^\Uk\i where Uk is 
the fc-th component of the Fourier spectrum of the displacement field ii„. Same parameters as in Fig. [5] 
Right panel; Velocity (triangles) of travelling multibreathers versus the wavenumber of the forcing. The 
solid line is the group velocity of linear waves. 


^ = 








7 = 


2 2 
— (jJZ 


h = 


■/! 7 f 

V 8(w2-w2)3/2-^’ 


and choosing a reference frame moving with velocity v = dujp/dp = sinp/ojp, Eq. (l4^ reduces 
to the well studied “externally” driven (or ac driven) damped NLS equation [66167) : 

d'P 92 ^ 

^ - he** . (50) 

Exact soliton solutions of this equation can be obtained for 7 ' = 0, see Eqs. (37-40) of Ref. [HH] , 
Moreover, multisoliton solutions are also derived in Ref. m- What we observe in Fig. [5] might 
well be a superposition of such solutions to form a train of “intrinsically localized” structures. 
However, one should bear in mind that NLS solutions can describe only low amplitude states. 
Therefore, they can be only a rough approximation of the pattern displayed in Fig. [SJ which 
shows high amplitude localized peaks. 

In the right panel of Fig. [ 6 l we plot the speed of the travelling multibreather as a function 
of the wavenumber of the forcing p-mode, which compares well with the group velocity of 
the corresponding linear waves, showing that nonlinear effects are negligible in this parameter 
range. 


3.3 Driving by one end: Ordinary and Bandgap Transmission 

To simulate the effect of an impinging wave, we impose to the /3-FPU chain [£ = 1 in Eq. (IT|)], 
the boundary condition 

uo{t) = Acosojt, (51) 

while free boundary conditions are enforced on the other side of the chain. 

In order to be able to observe a stationary state in the conducting regime, we need to 
steadily remove the energy injected in the lattice by the driving force. Thus, we damp a certain 
number of the rightmost sites (typically 10 % of the total) by adding a viscous term — 7 U„ 
to their equations of motion. A convenient indicator to look at is the averaged energy flux 
j = where the local flux is given by the following formula [55] 

jn = ^{Un+ Un+l) [un+1 “ Mn + {Un+1 “ ■ 


(52) 
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Time averages of this quantity are taken in order to characterize the insulating (zero flux)/con¬ 
ducting (non zero flux) state of the system. 


3.3.1 In-band driving: nonlinear phonons 

For illustration, we first discuss the case when the driving frequency is located inside the phonon 
band. Although trivial, this issue is of importance to better appreciate the fully nonlinear 
features described later on. 

Under the effect of the driving m, we can look for extended quasi-harmonic solutions 
(nonlinear phonons) of the form 

Un = Acos(kn — Lot). (53) 

We consider the semi-infinite chain, so that k varies continuously between 0 and 27r. The non¬ 
linear dispersion relation can be found in the rotating wave approximation (see e.g. Ref. m)- 
Neglecting higher-order harmonics, it reads 

LOQ{k, A) = 2(1 — cosk) -I- 3(1 — cos k)‘^A^. (54) 

Thus the nonlinear phonon frequencies range from 0 to the upper band-edge wo(7r, A) > 2. 

If we simply assume that only the resonating phonons whose wavenumbers satisfy the con¬ 
dition 

Lo = ujQ{k, A) (55) 

are excited, we can easily estimate the energy flux. Neglecting, for simplicity, the nonlinear 
force terms in the definition of the flux (IH^ . we have 

/ = iu(fc, A)a;^A^ (56) 

where v is the group velocity as derived from dispersion relation (1541) . This simple result is 
in very good agreement with simulations, at least for small enough amplitudes (see left panel 
of Fig. [21). For A > 0.15, the measured flux is larger than the estimate (l56l) . indicating that 
something more complicated occurs in the bulk (possibly, a multiphonon transmission) and 
that higher-order nonlinear terms must be taken into account. 




Fig. 7. Left panel: Average energy flux vs. driving amplitude for in-band forcing, co = 1.8, 7 = 5. Data 
have been averaged over 10® periods of the driving. The inset is an enlargment of the small-amplitude 
region and the dashed line is the single nonlinear phonon approximation (1561) . Right panel: Average 
energy flux vs. driving amplitude for out-band forcing, co = 3.5, 7 = 5. Data have been averaged over 
2 10® periods of the driving for a chain of A = 512 particles. 
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3.3.2 Out-band driving: supratransmission 


Let us now turn to the more interesting case in which the driving frequency lies outside the 
phonon band, oj > a;o(7r, 0) = 2. In a first series of numerical experiments, we have initialized 
the chain at rest and switched on the driving at time t = 0. To avoid the formation of sudden 
shocks [70], we have chosen to increase smoothly the amplitude from 0 to the constant value A 
at a constant rate, i.e. 


uq = Acos{ujt) 



-t/ri 


(57) 


where typically we set ti = 10. 

At variance with the case of in-band forcing, we observe a sharp increase of the flux at a 
given threshold amplitude of the driving, see right panel of Fig. [T] This phenomenon has been 
denoted as nonlinear supratransmission EH] to emphasize the role played by nonlinear localized 
excitations in triggering the energy flux. 

This situation should be compared with the one of in-band driving, shown in Fig. [T] where 
no threshold for conduction exists and the flux increases continuously from zero (more or 
less quadratically in the amplitude). Indeed, the main conclusion that can be drawn from the 
previous section is that there cannot be any amplitude threshold for energy transmission in 
the case of in-band forcing. Moreover, although at the upper band edge the flux vanishes, since 
it is proportional to the group velocity (see formula (1561) 1. it is straightforward to prove that 
it goes to zero with the square root of the distance to the band edge frequency. Hence, the 
sudden jump we observe in the out-band case cannot be explained by any sort of quasi-linear 
approximation. 

In the following, we investigate the physical origin of nonlinear supratransmission, distin¬ 
guishing the cases of small and large amplitudes. 

When the driving frequency is only slightly above the band (0<w — 2-Cl), one can resort 
to the continuum envelope approximation. Since we expect the zone-boundary mode = tt to 
play a major role, we let 

«« = (-l)"i[V'ne*"‘ + Ce-“‘]. (58) 

In the rotating wave approximation |69| and for slowly varying ijjn, one obtains from the FPU 
lattice equations the nonlinear Schrodinger equation (tpn —> 'ip{x,t)) [ZI] 


2ia;'0 = (w^ - A)tp - V’xx - 


(59) 




Fig. 8. Left panel: Snapshot of the local energy below the supratransmission threshold A = 0.15 < Ath 
for cj = 2.1 7 = 10. The initial condition is an envelope soliton dMl) with xo = -1-1.8. Right panel: 
Snapshot of particle displacements Un below the supratransmission threshold for a driving frequency 
w = 5.12 and a driving amplitude A = 0.5 < Ath = 2.05. One can observe, similarly to the left panel, 
that a moving discrete breather appears at the left boundary and propagates inside the bulk, leaving 
behind the static solution. 
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with the boundary condition 'ijj(0,t) = A. 

The well-known static single-soliton solution of Eq. (15^ corresponds to the family of enve¬ 
lope solitons (low-amplitude discrete breathers) 


= a(—1)" cos(ajt) sech — xo)a 


(60) 


with amplitude a 
condition to be 


y/( lj^ — 4)/6. The maximum of the soliton shape is fixed by the boundary 


xo = ± 


acosh(a/^) 

a '/6 


(61) 


In this approximation, we have two possible solutions: one with the maximum outside the chain, 
which is purely decaying inside the chain (minus sign in (1611) 1. and another with the maximum 
located within the chain (plus sign in (1611) 1. Overcoming the supratrasmission threshold corre¬ 
sponds to the disappearence of both solutions. Indeed, when the driving amplitude reaches the 
critical value Ath, given by 

La^ = 4 + 6Af^, (62) 


solution (1501) ceases to exist. 

We have investigated this issue by simulating the lattice dynamics with the initial conditions 
given by Eqs. dSOD and (|0T|). The evolution of the local energy En (see Eq. (|00)) 1 is shown in 
left panel of Eig. |8l The solution with the maximum inside the chain slowly moves towards the 
right and, eventually, leaves a localized boundary soliton (with maximum outside the chain) 
behind. The release of energy to the chain is non stationary and does not lead to a conducting 
state. 

The scenario drastically changes at the supratransmission amplitude Ath- The chain starts to 
conduct: a train of travelling envelope solitons is emitted from the left boundary (see left panel of 
Eig. [21). Here we should emphasize that the envelope soliton solution (1501) . which is characterized 
by the k = n carrier wave-number, has a zero group velocity. Thus, transmission cannot be 
realized by such envelope solitons. Instead, transmission starts when the driving frequency 
resonates with the frequency of the envelope soliton with carrier wave-number k = 7r(7V —2)/iV, 
next to the 7r-mode. However, as far as we consider a large number of oscillators {N = 500), 
we can still use expression (1621) for the 7r-mode frequency. 

The above envelope soliton solution (1601) is valid in the continuum envelope limit, and is 
therefore less and less accurate as its amplitude increases. Indeed, if the weakly nonlinear condi¬ 
tion is violated, the width of the envelope soliton becomes comparable with lattice spacing and, 
thus, one cannot use the continuum envelope approach. Eortunately, besides the slowly varying 
envelope soliton solution dsni, an analytic approximate expression exists for large amplitude 
static discrete breather solutions, which is obtained from an exact extended plane wave solution 
with “magic” wave-number 2tt/3 [72] 


Un 


a(—I)" cos [wB[a)t] cos 




(63) 


if |(7rn/3) ± xo| < 7r/2 and = 0 otherwise. 
Here xq is defined as follows 


xo = acos(7l/a), 


(64) 


where A is the driving amplitude. The breather frequency uisia) depends on amplitude a as 
follows 


— 1.03 


A/37r2(4 -|- Oa^) 
4/^(5) ’ 


(65) 


where K{s) is the complete elliptic integral of the first kind with argument s = 3a/y/2{9a'^ + 4) 
and the factor 1.03 takes into account a rescaling of the frequency of the “tailed” breather m 
(see also m)- As previously for the case of the envelope soliton solution, we perform a numerical 
experiment where we put initially on the lattice the breather solution of formula (I63|) . Choosing 









18 


Will be inserted by the editor 


the plus sign in this expression, we do not observe any significant transmission of energy inside 
the chain. Instead, the minus sign causes the appearance of a moving breather, which travels 
inside the chain leaving behind the static breather solution with plus sign. Right graph in Fig. [5] 
presents this numerical experiment. 

The static breather solution (1551) ceases to exist if the driving amplitude exceeds the thresh¬ 
old Ath given by the resonance condition 

uj = ujB{Ath)- (66) 

Above this threshold the supratransmission process begins via the emission of a train of moving 
breathers from the boundary, exactly as it happens in the case of small amplitudes. It should be 
mentioned again that the transmission regime is established due to moving discrete breathers. It 
has been remarked m that discrete breathers are characterized by quantized velocities, while 
their frequency is given by the same formula (16511 . This explains why one can use resonance 
condition dMl) for the static discrete breather solution (l63l) to define the supratransmission 
threshold in the large amplitude limit. 


3.3.3 Supratransmission threshold: numerical test 

To check these predictions, we have performed a numerical determination of Ath for several 
values of w, starting the chain at rest. This is accomplished by gradually increasing A and 
looking for the minimal value Ath for which a sizeable energy propagates into the bulk of the 
chain. At early time, the scenario is qualitatively similar to the one shown in the left graph 
of Fig. [9] Later on, the interaction of nonlinear and quasi-linear modes and their “scattering” 
with the dissipating right boundary establishes a steady energy flux into the chain. 

As seen in the right graph of Fig. [51 formulae (1551) [with definition (IHbll ] and (1551) (see the 
inset) are in excellent agreement with simulations for large A > 2 and small A < 1 amplitudes, 
respectively. The accuracy of the analytical estimate in formulae (1551) and (I55|) is of the order 
of few percents, at worst, in the intermediate amplitude range. 

For comparison, we have checked that the supratransmission threshold is definitely not 
associated with the quasi-harmonic waves with nonlinear dispersion relation (1541) . If this were 




Fig. 9. Left graph: Snapshot of the local energy at the transmission threshold A = 0.253 ~ ^th 
for oj = 2.1 7 = 10. The initial condition is the envelope soliton IMl) with xo ss 0. Right graph: 
Comparison between analytic estimates and numerical values of threshold amplitudes vs. the driving 
frequency. Main plot: the full dots are the numerical values of Ath and the solid line is a plot of formulae 
(1651661) . which are valid for large amplitudes. The inset shows an enlargement of the small Ath region, 
in order to illustrate the accuracy of the small-amplitude approximation (1621) (dotted line). 
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the case, the transmission should start when the oscillation amplitude reaches the value for 
which the resonance condition lo = uJo{k,A) holds. As u}o{k,A) is maximal for k = tt, we can 
get the expression for the threshold value from the relation w = Ath), he. 

a;2 = 4+12A?^. (67) 

The amplitude values one obtains from Eq. (1671) are far away from the numerical values and 
we don’t even show them in the right graph of Fig. [S] This is a further confirmation that 
supratransmission in the FPU model originates from direct discrete breather generation as it 
happens in the cases of discrete sine-Gordon and nonlinear Klein-Gordon lattices |55] . 


4 Conclusions 

In this paper, we have presented a detailed analysis of the zone-boundary mode modulational 
instability for the FPU lattice in both one and higher dimensions. Formulas for the critical 
amplitude have been derived analytically and compare very well with numerics for all system 
sizes. The study of the process which leads to the formation of chaotic breathers can be extended 
to two dimension; the physical picture is similar to the one-dimensional case. 

All results on modulational instability of zone-boundary modes can be straightforwardly 
extended to other initial modes and, correspondingly, instability rates can be derived. This 
has already been partially done in Ref. m and compares very well with the numerical results 
by Yoshimura m- This author has recently reanalyzed the problem m to determine the 
growth rates for generic nonlinearities in the high energy region, obtaining exact results based 
on Mathieu’s equation. 

For many-modes initial excitations, it has been remarked that instability thresholds depend 
on relative amplitudes and not only on the total energy [S] . Although this makes the study of the 
problem extremely involved, we believe that a detailed study of some selected group of modes, 
which play some special role in FPU dynamics, could be interesting. The method discussed in 
this paper could be adapted to treat this problem. Historically, the first study is in the paper 
by Bivins, Metropolis and Pasta himself m, where the authors tackle the problem by studying 
numerically the instabilities of coupled Mathieu’s equations. 

In one-dimensional studies, a connection between the average modulation instability rates 
and the Lyapunov exponents has been suggested [mm. Recently m. high frequency exact 
solutions have been used in the context of a differential geometric approach |63j to obtain 
accurate estimates of the largest Lyapunov exponent. Similar studies could be performed for 
the two-dimensional FPU lattice and the corresponding scaling laws with respect to energy 
density could be obtained. 

Finally, let us point out the generic nature of the results derived for the driven-damped Anti- 
FPU scenario. In this connection further developments towards nonlinear supratransmission 
and bistability effects in various physical systems, such as magnetic thin films m. Josephson 
junction arrays m, quantum Hall bilayers m optical directional couplers m and waveguide 
arrays [79180] should be expecially mentioned. 
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